Calculation of Gauss Quadrature Rules

Several algorithms are given and compared for computing Gauss quadrature rules. It is shown that given the three term recurrence relation for the orthogonal polynomials generated by the weight function, the quadrature rule may be generated by computing the eigenvalues and first component of the orthornormalized eigenvectors of a symmetric tridiagonal matrix. An algorithm is also presented for computing the three term recurrence relation from the moments of the weight function. | Introduction. Most numerical integration techniques consist of approximating the integrand by a polynomial in a region or regions and then integrating the polynomial exactly. Often a complicated integrand can be factored into a nonnegative "weight" function and another function better approximated by a polynomial, thus /b rb N g(t)dt = I o>(t)f(t)dt « 2>y/(iy) . u •'a j=l Hopefully, the quadrature rule {w¡, íyjyLi corresponding to the weight function co(i) is available in tabulated form, but more likely it is not. We present here two algorithms for generating the Gaussian quadrature rule defined by the weight function when : (a) the three term recurrence relation is known for the orthogonal polynomials generated by w(£), and (b) the moments of the weight function are known or can be calculated. In [6], Gautschi presents an algorithm for calculating Gauss quadrature rules when neither the recurrence relationship nor the moments are known. 1. Definitions and Preliminaries. Let w(x) ^ 0 be a fixed weight function defined on [a, b]. For o>(x), it is possible to define a sequence of polynomials po(x), pi(.x), ■ ■ ■ which are orthonormal with respect to w(x) and in which pn(x) is of exact degree n so that P I w{x)pm{x)pn{x)dx = 1 when m = n , (1.1) = 0 when m ?¿ n . The polynomial pn(x) = kn TFi=i (x — ti) ,kn> 0, has n real roots a < h < í2 < Received November 13, 1967, revised July 12, 1968. * The work of the first author was in part supported by the Office of Naval Research and the National Science Foundation; the work of the second author was in part supported by the Atomic Energy Commission. ** Present address: Hewlett-Packard Company, Palo Alto, California 94304. 221 222 GENE H. GOLUB AND JOHN H. WELSCH ■ • • < í„ < b. The roots of the orthogonal polynomials play an important role in Gaussian quadrature. Theorem. Letf(x) £ C2N[a, b], then fb N j-(2A0 /■.-. (2N)\h where Wi = fciV+l k r+i ! (n ,(t s dpN(t)\ \ n ptf+i(tj)pN'(tj) \ dt I «=,«./ j = l,2,---,N, Thus the Gauss quadrature rule is exact for all polynomials of degree ;£ 2N — 1. Proofs of the above statements and Theorem can be found in Davis and Rabinowitz [4, Chapter 2]. Several algorithms have been proposed for calculating {wj, tj)Nj-i ;cf. [10], [11] . In this note, we shall give effective numerical algorithms which are based on determining the eigenvalues and the first component of the eigenvectors of a symmetric tridiagonal matrix. 2. Generating the Gauss Rule. Any set of orthogonal polynomials, {pj(.x)}¥-i, satisfies a three term recurrence relationship : (-1) PÂ%) = (fl& + bj)pj-i(x) Cjpj-ïix), j= 1,2, ••-,#; p-i(x) = 0, p0(x) = l, with a, > 0, Cj > 0. The coefficients {a¡, b¡, c¡\ have been tabulated for a number of weight functions u>(x), cf. [8]. In Section 4 we shall give a simple method for generating {ay, bj, Cj] for any weight function. Following Wilf [12], we may identify (2.1) with the matrix equation